function Phi = basis_function_loc(V, x, y, basis_type, derivative_type)
Jacobi = [V(2,:)-V(1,:); V(3,:)-V(1,:)]';
J_det = det(Jacobi);
x_hat = (Jacobi(2,2)*(x-V(1,1)) - Jacobi(1,2)*(y-V(1,2)))/J_det;
y_hat = (-Jacobi(2,1)*(x-V(1,1)) + Jacobi(1,1)*(y-V(1,2)))/J_det;
if 4 == nargin || derivative_type == "x"
    Phi = basis_function_ref(x_hat, y_hat, basis_type);
else
    dNdx = basis_function_ref(x_hat, y_hat, basis_type, "dx");
    dNdy = basis_function_ref(x_hat, y_hat, basis_type, "dy");
    switch derivative_type
        case "dx"
            Phi = cellfun(@(x,y) (Jacobi(2,2)*x - Jacobi(2,1)*y)/J_det, dNdx, dNdy, "UniformOutput", false);
        case "dy"
            Phi = cellfun(@(x,y) (-Jacobi(1,2)*x + Jacobi(1,1)*y)/J_det, dNdx, dNdy, "UniformOutput", false);
    end
end
end

%% basis_function_ref
function N = basis_function_ref(x, y, basis_type, derivative_type)
switch basis_type
    case "P1"
        N = cell(3,1);
        if (3 == nargin)||("x" == derivative_type)
            N{1} = 1 - x - y;
            N{2} = x;
            N{3} = y;
        elseif "dx" == derivative_type
            N{1} = -ones(size(x));
            N{2} = ones(size(x));
            N{3} = zeros(size(x));
        elseif "dy" == derivative_type
            N{1} = -ones(size(x));
            N{2} = zeros(size(x));
            N{3} = ones(size(x));
        end
    case "P1b"
        N = cell(4,1);
        if (3 == nargin)||("x" == derivative_type)
            N{1} = 1 - x - y;
            N{2} = x;
            N{3} = y;
            N{4} = -x.*y.*(27*x + 27*y - 27);
        elseif "dx" == derivative_type
            N{1} = -ones(size(x));
            N{2} = ones(size(x));
            N{3} = zeros(size(x));
            N{4} = -27*y.*(2*x + y - 1);
        elseif "dy" == derivative_type
            N{1} = -ones(size(x));
            N{2} = zeros(size(x));
            N{3} = ones(size(x));
            N{4} = -27*x.*(x + 2*y - 1);
        end
    case "P2"
        N = cell(6,1);
        if (3 == nargin)||("x" == derivative_type)
            N{1} = 2*x.^2 + 4*x.*y - 3*x + 2*y.^2 - 3*y + 1;
            N{4} = -4*x.*(x + y - 1);
            N{2} = x.*(2*x - 1);
            N{5} = 4*x.*y;
            N{3} = y.*(2*y - 1);
            N{6} = -4*y.*(x + y - 1);
        elseif "dx" == derivative_type
            N{1} = 4*x + 4*y - 3;
            N{4} = 4 - 4*y - 8*x;
            N{2} = 4*x - 1;
            N{5} = 4*y;
            N{3} = zeros(size(x));
            N{6} = -4*y;
        elseif "dy" == derivative_type
            N{1} = 4*x + 4*y - 3;
            N{4} = -4*x;
            N{2} = zeros(size(x));
            N{5} = 4*x;
            N{3} = 4*y - 1;
            N{6} = 4 - 8*y - 4*x;
        end
    otherwise
        error("Invalid basis type.");
end
end